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This  paper  presents  a  fixed-time  glideslope  guidance  algorithm  that  is  capable  of 
guiding  the  spacecraft  approaching  a  target  vehicle  on  a  quasi-periodic  halo  orbit  in  real 
Earth-Moon  system.  To  guarantee  the  flight  time  is  fixed,  a  novel  strategy  for  designing 
the  parameters  of  the  algorithm  is  given.  Based  on  the  numerical  solution  of  the 
linearized  relative  dynamics  of  the  Restricted  Three-Body  Problem  (expressed  in  inertial 
coordinates  with  a  time-variant  nature),  the  proposed  algorithm  breaks  down  the 
whole  rendezvous  trajectory  into  several  arcs.  For  each  arc,  a  two-impulse  transfer  is 
employed  to  obtain  the  velocity  increment  (delta-v)  at  the  joint  between  arcs.  Here  we 
respect  the  fact  that  instantaneous  delta-v  cannot  be  implemented  by  any  real  engine, 
since  the  thrust  magnitude  is  always  finite.  To  diminish  its  effect  on  the  control,  a  thrust 
duration  as  well  as  a  thrust  direction  are  translated  from  the  delta-v  in  the  context  of  a 
constant  thrust  engine  (the  most  robust  type  in  real  applications).  Furthermore,  the 
ignition  and  cutoff  delays  of  the  thruster  are  considered  as  well.  With  this  high-fidelity 
thrust  model,  the  relative  state  is  then  propagated  to  the  next  arc  by  numerical 
integration  using  a  complete  Solar  System  model.  In  the  end,  final  corrective  control  is 
applied  to  insure  the  rendezvous  velocity  accuracy.  To  fully  validate  the  proposed 
guidance  algorithm,  Monte  Carlo  simulation  is  done  by  incorporating  the  navigational 
error  and  the  thrust  direction  error.  Results  show  that  our  algorithm  can  effectively 
maintain  control  over  the  time-fixed  rendezvous  transfer,  with  satisfactory  final 
position  and  velocity  accuracies  for  the  near-range  guided  phase. 

©  2012  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

The  near-earth  Rendezvous  and  Docking  (RVD)  pro¬ 
blem  has  been  studied  in  depth  [1-9]  since  the  1960s,  but 
not  too  much  progress  has  been  done  in  the  rendezvous 
problem  in  multi-body  realm.  However,  in  view  of  the 
great  potentials  of  LI  and  L2  in  Earth-Moon  system  (such 
as  the  portal  of  the  IPS  (Interplanetary  Superhighway) 
[10]  and  space-based  depot  for  lunar  exploration  [11]), 
the  problem  of  libration  point  orbit  rendezvous  should  be 
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well  addressed.  Unlike  the  near-earth  case,  practical 
orbits  around  libration  points  usually  have  very  large 
scales,  both  in  time  and  space.  One  of  the  most  frequently 
used  families  is  the  halo  orbit.  In  the  Circular  Restricted 
Three-Body  Problem  (CRTBP),  there  are  purely  periodic 
halo  orbits;  in  the  full  ephemeris  model,  only  quasi- 
periodic  halo  orbits  stand  [12],  which  we  will  refer  to  as 
real  halo  orbits  (or  halo  orbits  for  short)  in  this  paper.  The 
exact  determination  of  halo  orbits  has  to  resort  to  numer¬ 
ical  methods.  Moreover,  these  orbits  are  essentially  very 
unstable  due  to  the  high  nonlinearity  of  the  problem, 
posing  great  challenge  to  the  guidance  and  control  system 
design. 

There  are  already  some  works  focusing  on  the  halo-to- 
halo  transfers,  which  could  be  viewed  as  one  type  of  halo 
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orbit  rendezvous.  Hiday-Johnston  and  Howell  [13-15] 
applied  primer  vector  theory  [16]  to  the  transfer  design 
between  halo  orbits  in  Elliptical  Restricted  Three-Body 
Problem  (ERTBP),  and  found  some  local  cost-optimal 
results.  Gomez  et  al.  [17]  developed  two  methods  to  solve 
the  same  problem  in  CRTBP,  based  on  the  phase  space 
geometry  and  Floquet  theory,  respectively.  The  obtained 
results  were  compared  with  those  of  Hiday-Johnston  and 
Howell  [14],  Roberts  [18]  introduced  two  impulsive  stra¬ 
tegies  for  halo-to-halo  and  halo-to-Lissajous  transfers  for 
the  SOHO  (Solar  Heliospheric  Observatory)  mission.  Voile 
[19]  studied  the  halo  transfer  problem  in  the  Sun-Earth/ 
Moon  system  using  optimal  control  theory,  and  finite 
thrust  engine  with  variable  specific  impulse  was  assumed. 
In  general,  the  aforementioned  literature  takes  no  account 
of  either  the  gravitational  forces  other  than  from  the 
primaries,  or  errors  during  the  transfer  (such  as  the 
navigational  error  and  the  implementation  error,  or  control 
error),  thus  belongs  to  the  trajectory  design  or  mission 
planning  topics. 

The  scope  of  this  paper  falls  within  the  guidance  and 
control  problem  of  the  real  halo  orbit  rendezvous.  Only  a 
few  works  are  found  dealing  with  this  subject.  Jones  [20] 
has  done  some  pioneering  work  on  terminal  phase  ren¬ 
dezvous  in  CRTBP.  The  targeting  law  proposed  was  in  fact 
a  two-impulse  transfer  scheme,  which  is  only  applicable 
for  small  amplitude  halo  orbits  and  small  rendezvous 
ranges.  Marinescu  et  al.  [21,22]  studied  the  minimum 
propellant  optimal  low-thrust  rendezvous  problem,  but 
the  target  was  the  libration  point  itself,  and  only  planar 
results  were  given.  On  the  other  hand,  research  on  the 
near-earth  rendezvous  guidance  problem  is  abundant, 
among  which,  Hablani  et  al.  [23]  proposed  algorithms 
based  on  traditional  glideslope  guidance  with  astronauts 
in  the  loop,  for  spacecraft  to  approach,  to  fly  around,  and 
to  depart  from  a  target  vehicle  in  a  near-earth  circular 
orbit.  The  closed-form  solution  of  the  linear  Clohessy- 
Wiltshire  (C-W)  equations  [1]  was  utilized  and  the 
algorithms  were  based  on  the  impulsive  assumption. 
The  flight  time  in  Hablani’s  work  is  a  free  parameter  with 
no  special  attention  paid  to. 

In  view  of  the  importance  of  the  flight  time  in  RVD 
mission  (especially  with  man  involved),  the  purpose  of 
this  study  is  to  develop  an  algorithm  for  time-fixed 
missions  to  guide  the  spacecraft  to  approach  a  target 
vehicle  in  a  real  halo  orbit,  extending  the  near-earth  RVD 
to  multi-body  realm.  To  this  end,  several  reference  frames 
and  dynamical  models  are  introduced  in  Section  2. 
A  practical  constant  thrust  model  is  given  in  Section  3, 
having  considered  the  ignition  and  cutoff  delays  as  well  as 
a  linear  thrust  establishing  and  vanishing  process.  The 
guidance  algorithm  is  presented  in  Section  4,  in  which  a 
new  scheme  for  designing  the  control  parameters  for  a 
given  flight  time  is  derived.  In  Section  5,  impacts  of  these 
parameters  on  the  algorithm  performance  are  analyzed 
for  a  typical  standard  trajectory,  followed  by  Monte  Carlo 
simulations  in  consideration  of  navigation  and  control 
errors.  Results  have  reliably  justified  both  the  effective¬ 
ness  and  robustness  of  the  proposed  algorithm  in  real 
applications.  Appropriate  conclusions  are  drawn  in 
Section  6. 


2.  Equations  of  motion 

Spacecraft  are  usually  designated  chaser  and  target  in 
the  rendezvous  scenario,  in  which  the  target  is  intended 
to  follow  a  ballistic  trajectory  (such  as  a  halo  orbit  of 
interest)  with  infrequent  orbit  maintenance  maneuvers. 
Control  is  only  applied  to  the  chaser  to  approach  or  depart 
from  the  target. 

By  referring  to  the  near-earth  case  [7-9],  we  divide  the 
halo  orbit  rendezvous  into  three  consecutive  phases:  far- 
range  guided  phase,  near-range  guided  phase,  and  final 
approach  phase.  In  the  first  phase,  no  relative  navigation 
is  established  between  the  two  spacecraft,  and  the  chaser 
uses  absolute  navigational  means  to  approach  the  target. 
In  the  second  phase,  relative  measuring  information  is 
available  and  relative  motion  is  usually  considered  in  the 
guidance  methodologies.  Lastly,  the  final  approach  phase 
will  deliver  the  chaser  to  the  docking  position  through  a 
distance  of  hundreds  of  meters,  along  which  path  the  six- 
degree-of-freedom  (both  the  translational  and  rotational 
motion)  control  must  be  applied. 

In  this  study,  only  the  near-range  guided  phase  is 
considered.  In  the  first  place,  several  reference  frames 
and  dynamical  models  are  introduced  as  follows. 


2.1.  Geocentric  J2000.0  coordinates  frame  and  related 
dynamical  model 

The  origin  of  the  geocentric  J2000.0  frame,  designated 
by  Oe-XYZ,  is  at  the  barycenter  of  the  Earth.  The  X-axis 
points  at  the  mean  J2000.0  vernal  equinox,  and  the  X-Y 
plane  coincides  with  the  mean  equatorial  plane. 

The  equations  of  motion  described  in  Oe-XYZ  is  given 
by 

f  =  -^r+/M+/s+/P+^n  (1) 

where  r  =  §fif#  r  is  the  geocentric  position  vector,  /xE  the 
earth  gravitational  constant  (3.986004328969E14  m3 
s  2),  fM  the  lunar  gravitational  acceleration,  fs  the  solar 
gravitational  acceleration,  fP  the  sum  of  planetary  grav¬ 
itational  accelerations,  F  the  thrust  magnitude,  m  the 
spacecraft  mass,  and  n  the  thrust  direction  vector.  The 
accelerations  induced  by  the  celestial  bodies  are  calcu¬ 
lated  by  the  Newton’s  gravitational  law,  with  the  posi¬ 
tions  of  the  planets  given  by  the  JPL  DE405  ephemeris. 
Note  that  this  frame  is  not  inertial;  therefore,  the  impact 
of  the  motion  of  the  Earth  should  be  implicitly  removed 
from  those  gravitational  force  items. 

The  equation  related  to  mass  consumption  is  described 
in  [24] 

rh  =  -F/c  (2) 

where  c  is  the  exhaust  velocity. 

By  numerically  integrating  Eqs.  (1)  and  (2),  geocentric 
states  of  the  two  spacecraft  are  obtained,  from  which  the 
relative  states  are  further  calculated. 
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2.2.  Earth-Moon  LI  rotating  reference  of  frame 

Earth-Moon  LI  rotating  frame  (or  rotating  frame  for 
short),  designated  by  L-xry^r,  has  its  origin  at  the 
geometric  LI  point  which  is  located  between  the  two 
primaries.  The  xr-axis  is  along  the  vector  from  the  Earth  to 
the  Moon,  and  the  xryr  plane  is  the  instantaneous  lunar 
orbital  plane  with  zr-axis  as  its  normal  unit  vector. 

This  frame  is  established  by  using  the  real-time  posi¬ 
tions  of  the  Earth  and  the  Moon  read  from  JPL  DE405 
ephemeris  data.  Therefore,  all  the  irregular  motions  of 
both  the  Moon  and  the  Earth  are  implicitly  considered. 
The  quasi-halo  orbit  used  in  Section  5  is  expressed  in 
this  frame. 

2.3.  RVD  frame  and  relative  equations  of  motion 

The  RVD  frame,  designated  by  o—xyz,  is  chosen  to  have 
its  origin  at  the  barycenter  of  the  target,  and  the  three 
axes  identical  to  the  J2000.0  triad.  For  halo  orbit  rendez¬ 
vous,  much  convenience  with  the  coordinate  transforma¬ 
tion  is  brought  about  by  this  definition,  since  no  angular 
velocity  is  introduced. 

The  linearized  equations  of  relative  motion,  previously 
adopted  in  formation  flying  at  Lagrange  points  [25],  are 
employed  in  this  work  to  describe  the  dynamics  of  the 
chaser  relative  to  the  target  in  inertial  coordinates  of 
RTBP,  which  are 


where 


E(t)  =  —(ci  +c2)I3+3ci  e1TeJT+ 3c2e2TelT 
Cl  =  Hi  Ikir(t)  |  “3,C2  =  Hi  Ikzrft)  ||  “3 

and  /<,  is  the  gravitational  constant  of  the  larger  primary 
(3.986004328969E14  m3  s  2  for  the  Earth),  p2  the  gravita¬ 
tional  constant  of  the  smaller  primary  (4.902800582148E12 
m3  s  2  for  the  Moon),  r1T  the  target  position  relative  to 
larger  primary,  r2T  the  target  position  relative  to  smaller 
primary,  en  the  unit  vector  of  r1T,  e2J-  the  unit  vector  of 
r2T,  and  /3  a  3  x  3  identity  matrix. 

Note  that  the  time  variation  in  E  (or  A)  is  due  to  the 
relatively  slow  change  in  the  location  of  the  target  relative 
to  the  primaries.  Also  notice  that  the  relative  motion 
model  (3)  imposes  no  constraints  on  the  type  of  the 
spacecraft  orbit,  which  means  it  applies  to  all  the  orbit 
types  in  the  three  body  problem.  This  model  will  be 
utilized  in  the  proposed  guidance  algorithm  in  Section  4. 

3.  Constant  thrust  model 

Engines  of  constant  thrust  are  the  most  common  and 
robust  type  in  practical  applications.  Ideally,  it  has  a  step 
function  thrust  profile,  i.e.,  the  thrust  jumps  from  zero  to 
the  nominal  magnitude  in  an  instant  or  vice  versa. 
However,  the  actual  situation  is,  although  holding  con¬ 
stant  for  most  of  the  burn  time,  the  thrust  profile  has, 
inevitably,  both  ignition  and  cutoff  delays  and  a  contin¬ 
uous  rise-and-fall  process,  which  is  usually  difficult  to 


model  accurately  because  the  establishments  of  both  the 
propellants  pipe  flow  and  the  combustion  gas  filling  are 
highly  nonlinear.  Therefore,  a  constraint  of  minimum 
thrust  time  is  induced,  that  is,  a  null  or  insufficient 
magnitude  will  be  exerted  if  the  commanded  thrust  time 
is  too  transient,  which  will  surely  have  influences  on  the 
guidance  effect  (see  Section  5.1 ).  In  this  work,  a  simplified 
(but  practical  enough)  type  of  the  thrust  model  [26]  is 
used,  as  shown  in  Fig.  1. 

Relations  between  the  time-related  quantities  are 
given  by 

At,  =  (10/9)f(l0  +  tid  (4) 

At2  =  (10/9)tj0  +  tsd  (5) 

where  At,  is  thrust  build-up  time,  At2  the  thrust  tail-off 
time,  tg0  the  time  from  zero  to  90%  full  thrust,  tj0  the  time 
from  100%  to  10%  full  thrust,  tid  the  ignition  delay,  and  tsd 
the  cutoff  delay.  Note  that  At]  gives  the  minimum  time 
required  to  reach  the  nominal  thrust  magnitude. 

Fig.  1  shows  the  case  in  which  the  commanded  thrust 
time  (denoted  by  Atc)  is  sufficiently  long  and  the  nominal 
magnitude  can  be  achieved.  On  the  other  hand,  two  cases 
will  be  incurred  if  the  thrust  time  is  too  short,  distin¬ 
guished  by  comparing  with  tid  and  At, :  (1)  null  thrust,  if 
A tc<Ud,  (2)  insufficient  magnitude  other  than  zero,  if 
tid  <  Atc  <  At].  As  for  the  latter  case  (shown  in  Fig.  2),  we 
assume  the  thrust  profile  follows  the  same  decreasing 
rate  as  the  normal  one  shown  in  Fig.  1. 

In  conclusion,  the  thrust  in  the  ignition  and  the  cutoff 
phase  (denoted  by  F,  and  Fc),  respectively,  are  expressed 

control  | 


F/F 

1.0 

F,/Fm 


Fig.  2.  Thrust  profile  in  insufficient  thrust  time  case. 
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as  follows 

f  0  ?i  <  tid 

Fi  =  <  (0.9/fao)(r1-tid)Fm  tid  <  ti  <  Ati  (6) 

[Fm  Ti>Ati 

(  Fs  z2  <  tsd 

Fc  =  }  Fs-(0.9/t'w)(r2-tsd)Fm  tsd  <z2<  tr+tsd  (7) 

(o  r2>tr+tsd 

where  Ti=t~t;,T2=t-ts,  tr  =  (tj0/0.9)(Fs/Fm),  Fm  the 
nominal  thrust  magnitude,  and  Fs  the  thrust  when  cutoff 
order  is  given. 

4.  Constant  thrust  time-fixed  glideslope  algorithm 

Now  we  consider  the  time-fixed  rendezvous  problem. 
Assume  that  the  chaser  is  located  at  an  initial  position  r0 
with  an  initial  velocity  i/0  at  time  t0.  The  problem  is  how 
to  guide  it  to  r*  with  v*  within  a  given  flight  time  T,  where 
rf  is  the  chaser’s  desired  final  position  (could  be  a  desired 
parking  node  near  the  target  or  the  target  itself),  and  v* 
the  chaser’s  desired  final  velocity.  The  basic  idea  is  to 
break  down  the  trajectory  into  smaller  arcs,  which  are  to 
be  implemented  sequentially  using  the  two-impulse 
transfer  scheme  [20]  that,  though  highly  idealized  and 
unfit  for  a  real  flight,  provides  good  insight  for  the 
glideslope  guidance  algorithm.  Each  individual  arc  is 
associated  to  an  equally  spaced  guidance  period  (At), 
within  which  time  the  thruster  fires  along  a  fixed  direc¬ 
tion  for  a  certain  amount  of  time.  To  get  a  basic  idea,  see 
the  sketch  of  the  algorithm  depicted  in  Fig.  3. 

To  be  better  elaborated,  the  whole  algorithm  is  divided 
into  three  major  components:  the  design  of  the  time-fixed 
control  parameters,  the  glideslope  guidance  procedure, 
and  the  implementation  of  the  velocity  increments.  The 
first  component  deals  with  how  to  incorporate  the  flight 
time  into  the  control  law  design;  the  second  one  gives  the 
procedure  to  implement  the  rendezvous  transfer  using 
the  parameters  previously  designed;  the  last  one  focuses 
on  the  practical  thrusting  issues  when  finite  constant 
thrust  is  adopted. 


4.1.  Time-fixed  control  law  design 

Before  giving  the  procedure  of  the  algorithm,  we  first 
would  like  to  present  the  design  of  the  control  parameters 
that  are  to  be  used  in  the  algorithm.  Define  the  vector 
from  r0  to  rf  as  the  commanded  path  p,  then  we  require 
the  expected  location  at  time  t  equal  to 

r  =  r*t-pp  (8) 

where  p  =  p/||p||,  and  p  is  the  distance-to-go  along  p, 
whose  time  derivative  is  about  to  be  designed.  It  is 
demanded,  out  of  safety  consideration  in  RVD  mission, 
that  as  p  diminishes,  the  speed  p  must  diminish  with  it. 
Suppose  the  boundary  values  of  p0,  p0  and  pf  are  given 

Po  =  P(to)  =  if  r0i| 

P o  =  PV o)  <  0  (9) 

Pf  =  P(tf)  <  0 

where  p0  <  fij  according  to  the  safety  requirement,  and 
Po  is  the  rendezvous  distance.  Different  relationship 
between  p  and  p  can  be  postulated.  The  simplest  one  is 


the  linear  function,  given  by 

p  =  kp+b  (10) 

From  Eqs.  (9)  and  (10),  the  time  explicit  expression  for 
p  and  T  are  obtained  as 

p  =  (Po  +  b/k)exp(kt)-b/k  (11) 

T=ln(pf/p0)/k  (12) 

where 

k  =  Po~Pf/Po<0  03) 

b  =  Pf<  0  (14) 


Seen  from  Eqs.  (11)  and  (12),  both  p  and  T  are 
functions  of  p0  (a  known  parameter),  p0  and  pf,  which 
reasonably  makes  them  to  be  the  key  design  parameters 
of  our  interest. 

Unfortunately,  there  are  no  explicit  expressions  in 
which  p0(or  fif)  can  be  expressed  by  T,  which  means 
one  cannot  directly  solve  the  values  of  p0  and  pf  out  of  a 
given  T.  We  need  more  meaningful  parameters  to  bridge  a 
fixed  time  of  flight.  Note  that  the  total  number  of  arcs 
(denoted  by  N)  is  a  critical  parameter  in  the  sense  that  it  is 
associated  with  the  number  of  thruster  firings,  thus 
having  an  impact  on  both  the  energy  cost  and  guidance 
accuracy.  We  also  notice  that  the  p  value  of  the  final  arc 
(final  distance-to-go,  denoted  by  p*)  also  plays  an  impor¬ 
tant  role  in  the  control  performance.  By  introducing  these 
two  parameters,  both  with  clear  physical  meanings,  we 
present  an  algorithm  that  is  capable  of  solving  for  p0  and 
pf  with  a  given  rendezvous  distance  p0  and  the  flight  time 
T,  such  that  the  whole  guidance  law  can  be  established  in 
a  complete  way.  For  clarity,  we  put  it  in  the  form  of  a 
proposition,  and  give  the  deduction  in  the  appendix. 
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Proposition  1.  Assume  p0,  T,  N,  and  p*  ari 
and  Pf  are  to  be  solved.  Let 

s  given,  and  p0 

y  =  P*/Po 

(15) 

rl  =  Pf/Po 

(16) 

If  ye(0,l/N)  is  satisfied,  then  we  have 
p0  =  p0\n(ti)/[T(\-ri)] 

(17) 

Pf  =  rIPo 

(18) 

where  r\  is  the  solution  of 

^(N-i)/N  =  7+f7(1_7) 

(19) 

See  the  appendix  for  detailed  demonstration. 

Note  that  r\  is  a  fixed  unknown  in  Eq.  (19),  and  has  to 
be  solved  numerically  with  additional  efforts  (for 
instance,  by  Newton  method).  We  use 


MV-lf 

|N(i-y)J 


(20) 


as  its  initial  guess  (see  appendix  Eq.  (A5)),  which  enables 
a  fast  convergency  process. 

To  sum  up,  by  introducing  r],  variables  p0  and  pf  can 
now  be  explicitly  related  to  a  given  T.  Therefore,  inde¬ 
pendent  design  parameters  are  now  T,  N,  p0  and  p*,  based 
on  which  the  analysis  of  the  control  effects  will  be  done 
(see  Section  5). 


4.2.  Procedure  of  the  glideslope  guidance  algorithm 

Once  all  the  parameters  in  the  algorithm  have  been 
designed,  we  can  start  the  procedure  of  the  guidance 
process.  Although  for  brevity,  only  one  guidance  period  is 
discussed,  the  subscript  index  (i)  appearing  in  the  equa¬ 
tions  will  help  to  generalize  the  description  of  the  whole 
process. 

For  a  certain  arc  along  the  path,  the  i-th  (index  of  the 
thruster  firing,  i=l,...,N+l)  control  occurs  at  p(t,  , ) 
when  t,  !  =(i—  l)At  where  At=T/N  is  the  guidance  per¬ 
iod,  and  pushes  the  chaser  towards  p(q).  From  Eqs.  (8) 


and  (11),  the  corresponding  positions  are  given  by 
I’m  =  rf-[(p0+b/k)exp(kti-i)-b/k]p  (22) 

I)  =  r,_,  +  y  exp(kt,_i )p  (23) 

where  y  =  (p0+b/k)[\— exp(kAt)]  is  a  constant  value.  From 
Eq.  (3),  we  have 

r,  =  <P„  r,_,  +  0l2v+]  (24) 

vr  =  ^21*1-1  +  <f>22V+,  (25) 


where  vT1  represents  the  required  velocity  at  r,  ,  (navi¬ 
gational  position),  vr  the  arrival  velocity  at  r„  and 
d>mn,m,n=  1,2  the  four  3x3  submatrices  of  the  State 
Transition  Matrix  (STM)  <X>(t,-,ti  i),  which  is  the  solution 
of  the  following  matrix  Eq.  (26) 

{d>(t,_1,tM)  =  I6.  (26) 


Note  that  this  matrix  equation  has  to  be  numerically 
solved  with  the  dynamics  of  the  target  due  to  the  time 
variation  property  of  A. 

From  Eq.  (24),  we  have  the  required  velocity 
v+,  =  <f>72'(r,-<f>iiri-1)  (27) 

and  then  the  velocity  increment  at  rt_i 
Ay,- _v,.' ,-VC  j  (28) 

where  vr  1  is  the  chaser’s  arrival  velocity  given  by  the 
navigation  system  where  the  navigational  error  is 
introduced. 

Note  that  Ai/,  ,  cannot  be  implemented  by  a  finite 
constant  thrust  engine.  Suitable  transcription  of  this  velocity 
increment  is  required  (see  Section  4.3),  after  which,  the 
states  of  both  spacecraft  are  numerically  propagated  using 
the  full  Solar  System  ephemeris  until  the  next  guidance 
period  is  initialized. 

Special  attention  should  be  paid  to  the  final  corrective 
velocity  increment,  which  is  given  by 
Avn+i  =  v*-v))+1  (29) 

It  was  omitted  in  Ref.  [23],  but  based  on  the  results 
obtained  in  Sections  5.1  and  5.2,  we  consider  it  indispen¬ 
sable  for  diminishing  the  final  velocity  difference.  Note  that 
v))+1  is  not  the  navigational  velocity  but  the  one  predicted 
from  Eq.  (25).  This  is  because,  for  the  final  velocity  correc¬ 
tion,  the  thruster  should  work  before  reaching  the  final 
time,  and  only  if  the  value  of  ArN+1  is  given  can  the  thrust 
time  be  determined  (see  next  section).  Recall  that  p  is 
decreasing  as  time  goes  by,  therefore,  the  accuracy  of  the 
predicted  arrival  velocity  from  Eq.  (25)  is  growing,  making 
Eq.  (29)  a  good  approximation. 

4.3.  Implementation  of  velocity  increment 

In  the  context  of  the  finite  constant  thrust,  since  the 
magnitude  of  the  thrust  is  known  and  constant  (almost, 
such  as  the  thrust  model  we  use),  the  instantaneous  Avt  is 
transcribed  to  two  other  practical  variables:  the  thrust 
time,  and  the  direction. 

For  simplicity,  the  thrust  time  (denoted  by  At*.,  i  = 
1,...,N+1)  can  be  approximated  using  the  norm  of  Avt. 
Combining  Tsiollcovsky  equation  [24]  and  Eq.  (2),  we 
obtain  an  approximation  of  the  thrust  time 
At],  =  mc[l-exp(-AVj/c)]/Fm  (30) 

where  AVt  =  JAVj||  . 

Associated  with  the  thrust  time,  the  turn-on  time  of 
the  thruster  should  also  be  viewed.  Liu  [9]  used  a  1/2- 
time-ahead  scheme,  which  requires  the  prediction  of 
the  arrival  velocity.  In  this  study,  an  on-site  scheme  is 
used,  which  means  the  thruster  will  not  work  until 
the  thrust  time  is  given  by  the  guidance  system  based 
on  measurements  of  the  real-time  state.  In  the  last 
guidance  period,  however,  the  engine  should  be  ignited 
Atjy  ahead  in  order  to  correct  the  final  velocity  difference. 

Another  concern  is  the  thrust  direction.  We  simply 
hold  the  thrust  vector,  throughout  the  entire  thrust  time, 
to  be  aligned  with  the  direction  of  the  velocity  increment, 
i.e.,  tii=AVilAVi.  Errors  exist  in  holding  this  direction  (due 
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Required  velocity:  v*,  =  (r, 

N  I 

I  Velocity  increment:  Avf_,  =  v*_,  -  v~_,  L 


Translation  to  thrust  time  and  direction: 

A/j  =mc[l-exp(-A^/c)]/f.  n.  =  Av,/AK, 


State  propagation  using  JPL  ephemeris 


Fig.  4.  The  flowchart  of  the  glideslope  guidance  procedure. 

to  attitude  control  errors  and/or  thruster  installation 
errors)  and  should  be  taken  into  account  (see  Section  5.2). 

For  clarity,  we  give  a  flowchart  (see  Fig.  4)  of  the 
glideslope  guidance  procedure  as  a  conclusion  of  this  section. 


5.  Numerical  results 

Simulation  is  done  using  the  JPL  DE405  ephemeris, 
taking  into  account  all  the  gravitational  forces  from  Earth, 
Moon,  Sun  and  other  planets.  The  epoch  is  arbitrarily 
selected  to  be  56279.32388311  MJD  (18  Dec  2012 
07:46:23.500  UTCG).  A  location  on  a  quasi-periodic  halo 
orbit  (8000  km  of  z-amplitude,  in  the  real  Earth-Moon 
system)  is  employed  as  the  initial  target’s  position,  as 
shown  in  Fig.  5.  The  chaser’s  properties  are  listed  in 
Table  1,  where  all  the  thrust-related  values  are  normal 
to  the  current  engine  technology.  Without  losing  general¬ 
ity,  the  initial  and  final  relative  states  of  the  chaser,  not 
necessarily  on  a  halo  orbit,  are  given  in  Table  2,  where  D 
(D=pa  since  y=z=0)  is  the  rendezvous  distance  (the 
tested  range  is  100-2000  km  in  this  work,  coherent  to 
the  one  covered  by  the  typical  near-range  guided  phase 
in  RVD). 

Using  the  proposed  algorithm,  any  parameter  tetrad 
(T,N,p*,D)  would  give  a  rendezvous  orbit.  To  assess  the 
control  performance,  four  indices  are  surveyed,  which  are 
the  final  position  error,  the  final  velocity  error,  total  velocity 
increments,  and  mass  consumption,  defined  as  follows, 


mo  (kg) 
Fm(N) 
c  (m/s) 


3000 

1000 

2800 

0.002 

0.002 

0.010 


vx  (m/s) 
vy  (m/s) 
vz  (m/s) 


respectively 

ll^/IMIr/— */!  (3D 

lAv/l  =  llv/-vfll  (32) 

AV  =  J2Avi  (33) 

Am  =  rrif—mo  (34) 

where  V/  is  the  actual  final  velocity  of  the  chaser,  m0  the 
initial  mass,  rrif  the  final  mass. 


5.1.  Standard  trajectory  case 

In  the  standard  trajectory  case,  no  error  is  considered. 
The  objective  is  to  validate  the  feasibility  of  the  proposed 
guidance  algorithm. 
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First,  we  investigate  the  impact  of  D  and  T,  which  are 
closely  related  to  a  specified  mission.  As  a  matter  of  fact, 
the  shorter  the  flight  time  (usually  the  more  preferable  in 
RVD)  is,  the  more  difficult  it  will  be  for  the  rendezvous 
control.  In  view  of  the  possible  time-rigorous  missions, 
such  as  the  crew  retreat  from  a  space  station  in  emer¬ 
gency,  the  flight  time  we  consider  in  this  paper  will  be 
between  3-12  h  (flight  time  less  than  3  h  tends  to  give 
very  unfavorable  performance  for  most  cases  in  our  tests). 
Arbitrarily  set  N=5  and  p*=  50  m.  The  surfaces  in  Fig.  6 
show  the  control  performance  with  different  D—T  com¬ 
binations,  whereon  the  points  mark  the  calculated  cases. 
As  can  be  seen,  for  a  specific  D  (in  km),  different  T  will 
result  in  distinct  final  position  errors  ( ||  Aty  ||  values),  some 
of  which  even  reach  the  magnitude  of  10  km,  while  some 
are  less  than  1  m.  Very  large  ||  Aiy  ||  occurs  only  when  T  is 
comparatively  short,  which  means  an  over-rigid  require¬ 
ment  of  flight  time  will  exceed  the  thruster’s  ability.  As 
expected,  the  best  choice  of  T,  which  produces  the 
smallest  ||  Ary  || ,  grows  as  D  increases.  A  similar  situation 
happens  to  ||Aiy||,  and  in  all  cases  it  lays  below  4  mm/s, 
which  is  more  than  sufficient  for  a  real  mission.  As  for  AV 
and  Am,  they  both  decreases  smoothly  as  T  increases 
within  its  range.  This  fact  indicates  that  the  control 
accuracy  and  the  energy  cost  cannot  be  satisfied  most  at 
the  same  time.  A  compromise  between  guidance  accuracy 
and  energy  cost  has  to  be  made  in  determining  a  proper  T 
for  a  given  D. 

Note  that  three  unexpected  points  appear  in  the  upper 
two  subplots  in  Fig.  6,  whose  ||Aiy|j  values  have  shown 
abrupt  rises.  The  one  with  D=  1200  km  and  T=10h  is 
sorted  out  for  a  close  inspection.  Fig.  7  shows  the  detail  of 
this  instance,  in  which  the  orbit  in  RVD  frame  is  plotted  in 


Fig.  7a  and  b,  the  orbit  in  rotating  frame  in  Fig.  7c,  and  the 
thrust  and  mass  profiles  in  Fig.  7d.  Thicker  lines  indicate 
thrusting  arcs  (Fig.  7a-c),  while  the  circle  and  the  points 
(Fig.  7a  and  b)  signify  r*  and  the  expected  locations  along 
p,  respectively.  One  can  see  in  Fig.  7b,  the  last  but  one 
control  (where  there  is  an  inflection)  is  unable  to  draw  the 
orbit  towards  the  destination,  thus  a  relatively  large 
|| Ary  ||  has  been  produced.  Checking  the  thrust  history, 
we  know  the  reason  lies  in  that  the  corresponding  burn 
time  (AtN_i  =  0.0112  s)  is  too  short  to  reach  the  nominal 
thrust  magnitude,  as  shown  in  Fig.  7d  (at  t=8  h).  Similar 
occasions  happened  to  the  other  two  cases.  Therefore,  we 
can  conclude  that  the  insufficient  implementation  of  Alt, 
will  largely  hazard  the  position  accuracy,  and  should  be 
paid  special  attention  to  in  the  finite  thrust  guidance 
applications. 

Next,  consider  the  influence  of  N-T  combinations  on 
the  control  performance.  Calculations  are  done  with 
respect  to  typical  values  of  D=500km  and  p*=50m. 
Note  that  N  is  discrete  integers.  Results  are  given  in  Fig.  8. 
In  all  the  four  subplots,  curves  tend  to  gather  as  N 
increases,  implying  too  many  thruster  firings  are  unne¬ 
cessary  for  achieving  certain  performance.  ||Aiy||  varies 
greatly  in  magnitude,  and  its  minimum  happens  at  N=  5, 
T=5  h.  Though  varying  abruptly  for  some  N  as  well,  ||  Aiy  || 
is  kept  under  good  control  (all  below  1  mm/s,  and  curves 
for  N>  6  are  hardly  discernable).  In  contrast,  AV  and  Am 
change  much  regularly  with  N  and  T.  They  decrease  as  T 
increases,  and  behave  reversely  with  respect  to  N,  which 
means  too  short  a  flight  time  or  too  many  thruster  firings 
will  induce  larger  energy  cost. 

Last,  check  the  impact  of  p*  selections.  Other  three 
parameters  are  chosen  to  be  D= 500  km,  T=5  h,  and  N=  5. 


Fig.  6.  Control  effects  with  different  D-T  combinations  (N=5  and  p*=50  m). 
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Results  are  given  in  Fig.  9,  where  points  on  the  lines 
represent  the  calculated  instances.  One  can  see  that  ||  Af/|| 
shows  a  rising  trend  as  p*  grows,  with  the  exception  at 


p*=\ 00  m  where,  as  expected,  an  insufficient  delta-v 
implementation  occurs.  AV  and  Am  both  decrease  as  p* 
increases,  and  around  100  m/s  of  delta-v  difference  exists 


y.  Lian  et  al./Acta  Astronautica  79  (2012)  241-252 


249 


Fig.  9.  Impact  of  p*  selections  (D=500  km,  T=5  h,  at 


for  a  5  km  range  of  p*  which  indicates  that  p*  has  a  great 
impact  on  the  energy  cost. 

To  sum  up  this  subsection,  a  large  number  of  tests, 
which  are  grouped  by  the  independent  design  para¬ 
meters,  have  validated  the  proposed  algorithm’s  feasibil¬ 
ity.  Each  parameter  in  the  set  (T,  N,  p*  D)  plays  a  different 
role  in  affecting  the  overall  control  performance.  Although 
optimization  is  not  the  subject  of  this  work,  we  can 
still  conclude  that  Multi-Objective  Optimization  (MOO) 
method  should  be  adopted  to  get  the  best  combination  (in 
fact,  there  will  be  many  since  the  Pareto  frontier  is  a 
solution  set)  of  these  parameters  when  both  the  guidance 
accuracy  and  energy  cost  are  considered. 

5.2.  Monte  Carlo  analysis 

This  subsection  deals  with  the  robustness  of  the 
proposed  algorithm.  In  real  applications,  there  are  various 
errors  existing  and  the  orbit  will,  inevitably,  deviate  from 
its  standard  path  (the  designed  one).  We  are  interested  in 
the  performance  of  the  algorithm  when  errors  are  intro¬ 
duced.  Navigation  error  and  control  error  are  the  main 
concerns  in  this  study. 

Navigation  error,  including  the  position  error  and  the 
velocity  error,  comes  from  the  navigational  equipment 
and  methods,  which  is  assumed  to  be  in  regard  with  the 
distance  between  chaser  and  rj  (denoted  by  R,  varying 
with  time)  in  this  analysis.  In  addition,  thrust  direction 
error,  resulting  from  attitude  control  error  or  thruster 
mounting  error,  is  also  considered.  Both  errors  are 


Standard  deviation 

R<  2km 

R>  2  km 

<tPi  (m) 

0.1 

1 

<7vj  (m/s) 

0.001 

0.01 

<TAVj 

0.5% 

assumed  to  be  normally  distributed  with  zero  means, 
and  the  standard  deviations  are  given  in  Table  3,  where 
aPj  is  the  standard  deviation  of  navigational  position 
error,  aVj  the  standard  deviation  of  navigational  velocity 
error,  oaV/  the  standard  deviation  of  thrust  direction  error, 
and  j=x,yz. 

Monte  Carlo  analyses  based  on  300  trials  are  con¬ 
ducted,  with  respect  to  three  cases  that  consider  only  the 
navigational  error,  only  the  thrust  direction  error,  and  the 
both,  respectively.  The  baseline  standard  trajectory  is 
taken  from  Fig.  6,  with  D=500km,  7=5  h,  N=5,  and 
p*=50m.  Results  are  given  in  Table  4  and  Fig.  10.  In 
Table  4,  statistics  of  ||  Aiy  ||,  |AVfj|(  AV  and  Am  in  the  three 
cases  are  listed:  the  mean,  the  standard  deviation,  the 
maximum,  and  the  minimum  of  the  300  trials.  Values  of 
the  standard  orbit  are  also  given.  Fig.  10  shows  the 
statistical  distribution  of  the  third  case  in  which  both 
errors  are  considered. 

It  is  apparent  in  Table  4  that  the  navigation  error 
accounts  for  most  part  of  the  deviations.  |Aiy|  is  largely 
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Table  4 

Monte  Carlo  simulation  results  (300  trials). 


Case 

Statistics 

IIAiy!  Cm) 

IIAv/ll  (m/s) 

AV  (m/s) 

Am  (kg) 

Standard  orbit 

0.141 

8.088e— 4 

246.707 

253.016 

Navigational  error  only 

Average 

5.745 

0.002 

246.710 

253.019 

Std.  deviation 

2.507 

7.365e  — 4 

0.036 

0.035 

Maximum 

13.041 

0.004 

246.814 

253.121 

Minimum 

0.537 

2.616e— 4 

246.619 

252.929 

Average 

0.163 

8.090e— 4 

246.711 

253.020 

Std.  deviation 

0.040 

1.005e— 6 

0.003 

0.003 

Maximum 

0.455 

8.181e— 4 

246.719 

253.028 

Minimum 

0.133 

8.065e— 4 

246.701 

253.010 

5.703 

0.002 

246.709 

253.018 

Std.  deviation 

2.604 

7.916e— 4 

0.035 

0.035 

Maximum 

14.258 

0.005 

246.817 

253.124 

Minimum 

0.564 

6.244e  —  5 

246.619 

252.930 

AY  (m/s)  Am  (kg) 

Fig.  10.  Statistical  distribution  of  control  effect  indices  with  both  errors  considered. 


affected  by  the  introduced  errors,  with  the  maximum 
value  near  15  m,  while  the  minimum  around  0.6  m  when 
both  errors  are  considered.  In  contrast,  ||  Av/ 1| ,  AY,  and  Am 
are  only  slightly  influenced  by  the  errors.  ||Af/j|  is  always 
under  5  mm/s,  and  the  maximum  deviations  of  AV  and 
Am  are  also  less  than  0.1  m/s  and  0.1  kg,  respectively. 

Therefore,  it  can  be  concluded  that  the  last  firing  of  the 
algorithm  has  very  satisfactory  corrective  effect  on  the 
final  relative  velocity,  but  not  likely  on  the  position. 
In  contrast,  the  control  prior  to  the  last  one  is  very 
essential  for  the  final  position  accuracy.  We  should  also 
notice  that,  although  decreased  as  compared  with  the 
standard  case,  the  guidance  accuracy  is  still  sufficient  for 


initializing  a  final  approach  phase  that  will  finally  guide 
the  chaser  to  dock  with  the  target. 

6.  Conclusions 

The  proposed  constant-thrust  glideslope  algorithm 
provides  an  effective  approach  to  guide  the  spacecraft  to 
rendezvous  with  a  target  in  multi-body  realm.  It  is 
remarkable  that  a  new  strategy  for  designing  parameters 
with  the  fixed-time  constraint  is  presented,  which  will 
benefit  missions  that  are  time-sensitive.  Simulations 
indicate  that  the  final  position  accuracy  should  be  the 
primary  concern,  which  tends  to  be  easily  affected  by 
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insufficient  implementation  of  the  delta-v.  With  common 
errors  in  consideration,  the  method  can  reach  an  average 
control  accuracy  of  less  than  6  m  and  2  mm/s  (for  a 
rendezvous  range  of  500  km)  for  the  final  position  and 
velocity,  respectively,  making  it  a  good  candidate  for  the 
near-range  guided  phase  in  halo  orbit  rendezvous. 

Moreover,  the  core  of  the  method  is  quite  general.  It 
applies  not  only  to  the  approaching  operation  in  a 
rendezvous  mission,  but  retreating  or  fly-around  move¬ 
ments  as  well.  The  differences  lie  in  the  characteristics  of 
the  target  only.  Moreover,  the  dynamics  involved  in  the 
algorithm  add  no  constraints  on  either  the  type  of  the 
orbits  or  the  scope  of  the  three  body  system.  For  instance, 
it  can  be  employed  to  implement  the  Lissajous  orbit 
rendezvous  in  the  Sun-Earth/Moon  system. 


and  if  y-*l /N,  g(y)->  -0.  In  addition,  for  Vye(0,l/N), 
g'00  =  (l-Ny)(l-y)N-2  >  0 

Now  we  can  safely  draw  the  conclusion  that  Eq.  (19) 
will  always  have  only  one  solution  if  ye(0,l/N)  is 
satisfied. 

From  Eqs.  (12),  (15),  and  (16),  we  have 
T  =  p0  ln()7)/[/30(l  -»?)]  (A9) 

Then,  from  Eqs.  (A9)  and  (16),  p0  and  pj  are  given  by 

p0=Po 

Pf  =  rIPo 

which  are  Eqs.  (17)  and  (18),  respectively. 
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Appendix  A.  Demonstration  of  Proposition  1 


From  the  definition  of  p*.  we  have 
p*  =  p(tN_i)  (Al) 

Using  Eqs.  (11)  and  (Al), 

(p0+b/k)exp[(N-f  )kT /N]-b/k  =  p*  (A2) 

Substituting  Eqs.  (12),  (15)  and  (16)  into  Eq.  (A2)  yields 

^(N-D/N  =y+ >7(1-7) 

which  is  Eq.  (19),  where  17,7  e( 0,1)  by  their  definitions. 
Now  analyze  the  solutions  of  Eq.  (19).  Set 

fW  =  J7<N-1)/JV— »7(l  —y)—y  (A3) 

where  r/,y  e  (0,1 ).  Note  that  if  r\  -*  0/-»  —  y  <  0,  and  if  17  ->  1 , 
/-> 0,  thus  one  solution  of  Eq.  (A3)  can  be  insured  only  if 
f(i\)  >  0  (A4) 


where  fj  gives  the  maximum  of  f[})). 

Differentiating  Eq.  (A3)  and  setting /(t/)=0,  we  have 


■  \JtUT 

|_N(1  — y)J 


(A5) 


Recall  that  fj  e  (0,1).  Thus,  a  strengthened  constraint 
on  y  is  solved  as 

0  <  y  <  1  /N  (A6) 

Substituting  Eq.  (A5)  into  Eq.  (A4)  yields 

<(N--1)n-'/Nn  (A7) 

Now  check  the  inequality  of  Eq.  (A7).  Set 
g(y)  =  7(1  -7)N_1  -(N-  1)N_1  /Nn  (A8) 

and  g(y)  <  0  is  required.  From  Eqs.  (A6)  and  (A8),  we  see 
that  if  y->0, 

g(y)->-(N-lf~1/NN<0 
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